pacman::p_load(sf, raster, spatstat, sparr, tmap, tidyverse, stpp, skimr, gganimate, ggplot2, plotly, flexdashboard, shiny, DT,gridExtra, stars)Take_Home Exercise 1: Geospatial Analytics for Social Good: Application of Spatial and Spatio-temporal Point Patterns Analysis to discover the geographical distribution of Armed Conflict in Myanmar
1.0 Introduction
The study of armed conflicts in Myanmar has gained critical importance in understanding the geographical distribution and intensity of violence across different regions. Myanmar’s complex ethnic composition and ongoing civil strife make it a unique case for geospatial analysis. This project aims to apply spatial and spatio-temporal point pattern analysis methods to uncover the patterns of armed conflict between January 2021 and June 2024.
By leveraging conflict data from the Armed Conflict Location & Event Data (ACLED) and geospatial tools, we will focus on visualizing and interpreting conflict density through heat maps, Kernel Density Estimation (KDE), and advanced spatio-temporal analysis.
Our analysis will focus on four types of conflict events:
- Battles,
- Explosion/Remote violence,
- Strategic developments,
- Violence against civilians,
with particular attention paid to quarterly patterns in conflict occurrence.
2.0 Setting Up The Environment & Dataset
2.1 Installing the required Packages
Key Packages Used in the Project:
sf: Handles simple features in R, allowing for spatial data manipulation and analysis. It is crucial for reading and managing geospatial data like shapefiles (e.g., Myanmar’s administrative boundaries).raster: Used for raster-based spatial data manipulation, especially for working with raster maps, such as Kernel Density Estimation (KDE) results.spatstat: A powerful package for spatial point pattern analysis. It helps to analyze and visualize spatial point data, particularly for identifying clusters or patterns in armed conflict events.sparr: Builds onspatstatand focuses on performing spatial and spatio-temporal kernel smoothing, which will be crucial for KDE and heatmap creation.tmap: A thematic mapping package that will allow us to create maps, including KDE visualizations on an OpenStreetMap base.tidyverse: A collection of data manipulation packages likedplyr,ggplot2, andpurrr. It’s essential for data cleaning, manipulation, and visualization tasks.stpp: Used for spatio-temporal point pattern analysis, crucial for analyzing how conflict events evolve in both space and time.skimr: A quick and comprehensive tool to provide summaries and descriptive statistics for datasets, helping in the initial exploration.gganimate: Extendsggplot2to create animated visualizations. We can use this for animated time-series or evolving conflict maps.ggplot2: The core plotting package in R, essential for creating visualizations like time series plots and KDE heatmaps.plotly: Useful for creating interactive visualizations, allowing users to explore spatial data interactively (e.g., hover over points to see conflict details).pacman: is a package management tool in R designed to streamline the process of loading and installing packages.
2.2 Data-set involved in this topic
For this analysis, we use two key datasets:
2.2.1 ACLED Armed Conflict Data
Location & Event Data (ACLED)platform, which maintains an extensive record of conflict events globally. For this specific analysis, we will limit the dataset by filtering based on the following parameters to streamline data preparation and minimize the need for extensive data cleaning:
| Data Parameter | Filter Category |
|---|---|
| Date Range | From 01/01/2021 to 30/06/2024. |
| Event Type | 1. Battles 2. Violence Against Civilians 3. Explosions/Remote Violence 4. Strategic Developments |
| Country | Myanmar |
ACLED Configuration Image
Code to Import ACLED Dataset
ACLEDData <- read_csv("data/raw/aspatial/2021-01-01-2024-06-30-Myanmar.csv")Rows: 42608 Columns: 35
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (20): event_id_cnty, event_date, disorder_type, event_type, sub_event_ty...
dbl (15): year, time_precision, inter1, inter2, interaction, iso, latitude, ...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.2.2.1.1 Understanding the data set fields.
Referencing this ACLED Official codebook, this is the dataset that we are working with, not to bore you with the details are mainly interested in the following fields,
Event ID: Unique identifier for each conflict event.
Event Date: Date of occurrence.
Event Type: Type of conflict event (e.g., Battles, Remote Violence).
Latitude/Longitude: Coordinates of the event.
Fatalities: Number of fatalities resulting from the event.
Actors: The groups or individuals involved in the conflict (e.g., state actors, ethnic armed groups).
Admin Levels: Administrative region, district, and township where the event took place.
If you’re interested in the data set fields to explore more, here’s the full fields!
ACLED Full Table Fields Summary
| Fields name | Fields Description | Values |
| event_id_cnty | A unique alphanumeric event identifier by number and country acronym. This identifier remains constant even when the event details are updated. | E.g. ETH9766 |
| event_date | The date on which the event took place. Recorded as Year-Month-Day. | E.g. 2023-02-16 |
| year | The year in which the event took place. | E.g. 2018 |
| time_precision | A numeric code between 1 and 3 indicating the level of precision of the date recorded for the event. The higher the number, the lower the precision. | 1, 2, or 3; with 1 being the most precise. |
| disorder_type | The disorder category an event belongs to. | Political violence, Demonstrations, or Strategic developments. |
| event_type | The type of event; further specifies the nature of the event. | E.g. BattlesFor the full list of ACLED event types, see the ACLED Event Types table. |
| sub_event_type | A subcategory of the event type. | E.g. Armed clashFor the full list of ACLED sub-event types, see the ACLED Event Types table. |
| actor1 | One of two main actors involved in the event (does not necessarily indicate the aggressor). | E.g. Rioters (Papua New Guinea) |
| assoc_actor_1 | Actor(s) involved in the event alongside ‘Actor 1’ or actor designations that further identify ‘Actor 1’. | E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank. |
| inter1 | A numeric code between 0 and 8 indicating the type of ‘Actor 1’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). | 1, 2, 3, 4, 5, 6, 7, or 8. |
| actor2 | One of two main actors involved in the event (does not necessarily indicate the target or victim). | E.g. Civilians (Kenya)Can be blank. |
| assoc_actor_2 | Actor(s) involved in the event alongside ‘Actor 2’ or actor designation further identifying ‘Actor 2’. | E.g. Labor Group (Spain); Women (Spain)Can have multiple actors separated by a semicolon, or can be blank. |
| inter2 | A numeric code between 0 to 8 indicating the type of ‘Actor 2’ (for more, see the section Actor Names, Types, and ‘Inter’ Codes). | 0, 1, 2, 3, 4, 5, 6, 7, or 8. |
| interaction | A two-digit numeric code (combination of ‘Inter 1’ and ‘Inter 2’) indicating the two actor types interacting in the event (for more, see the section Actor Names, Types, and ‘Inter’ Codes). | E.g.3, 58 |
| civilian_targeting | This column indicates whether the event involved civilian targeting. | Either ‘Civilians targeted’ or blank. |
| iso | A unique three-digit numeric code assigned to each country or territory according to ISO 3166. | E.g. 231 for Ethiopia |
| region | The region of the world where the event took place. | E.g. Eastern Africa |
| country | The country or territory in which the event took place. | E.g. Ethiopia |
| admin1 | The largest sub-national administrative region in which the event took place. | E.g. Oromia |
| admin2 | The second largest sub-national administrative region in which the event took place. | E.g. ArsiCan be blank. |
| admin3 | The third largest sub-national administrative region in which the event took place. | E.g. MertiCan be blank. |
| location | The name of the location at which the event took place. | E.g. Abomsa |
| latitude | The latitude of the location in four decimal degrees notation (EPSG:32647). | E.g. 8.5907 |
| longitude | The longitude of the location in four decimal degrees notation (EPSG:32647). | E.g. 39.8588 |
| geo_precision | A numeric code between 1 and 3 indicating the level of certainty of the location recorded for the event. The higher the number, the lower the precision. | 1, 2, or 3; with 1 being the most precise. |
| source | The sources used to record the event. Separated by a semicolon. | E.g. Ansar Allah; Yemen Data Project |
| source_ scale | An indication of the geographic closeness of the used sources to the event (for more, see the section Source Scale). | E.g. Local partner-National |
| notes | A short description of the event. | E.g. On 16 February 2023, OLF-Shane abducted an unidentified number of civilians after stopping a vehicle in an area near Abomsa (Merti, Arsi, Oromia). The abductees were traveling from Adama to Abomsa, Arsi. |
| fatalities | The number of reported fatalities arising from an event. When there are conflicting reports, the most conservative estimate is recorded. | E.g. 3No information on fatalities is recorded as 0 reported fatalities. |
| tags | Additional structured information about the event. Separated by a semicolon. | E.g. women targeted: politicians; sexual violence |
| timestamp | An automatically generated Unix timestamp that represents the exact date and time an event was uploaded to the ACLED API. | E.g. 1676909320 |
2.2.2 Myanmar Administrative Boundaries (Shapefiles):
Obtained through Geonode Mimu, this shapefile helps us to build the map and set the boundary zone of each district of myanmar. This dataset provides the geographical boundaries of Myanmar’s administrative divisions, from the national level down to the township level. It is essential for mapping conflict events to specific regions.
Note
Myanmar has State, District and Township level, why District Level?
Choosing the district level over the township level for conflict analysis provides a better balance between detail and clarity. The district level allows us to capture broader regional trends without overwhelming the analysis with too many granular data points, as township-level data can be overly detailed. It also improves computational efficiency and makes visualizations clearer, while still offering enough specificity to reveal conflict hotspots. Additionally, population and auxiliary data are more readily available at the district level, making the analysis more consistent and manageable.
Code to Import Shapefile Dataset
M_State_Sf <- st_read(dsn="data/raw/geospatial/stateLevel", layer = "mmr_polbnda_adm1_250k_mimu_1") Reading layer `mmr_polbnda_adm1_250k_mimu_1' from data source
`C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial\stateLevel'
using driver `ESRI Shapefile'
Simple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84M_State_SfSimple feature collection with 15 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84
First 10 features:
OBJECTID ST ST_PCODE ST_RG ST_MMR PCode_V
1 1 Ayeyarwady MMR017 Region ဧရာဝတီတိုင်းဒေသကြီး 9.4
2 2 Bago MMR111 Region ပဲခူးတိုင်းဒေသကြီး 9.4
3 4 Chin MMR004 State ချင်းပြည်နယ် 9.4
4 5 Kachin MMR001 State ကချင်ပြည်နယ် 9.4
5 6 Kayah MMR002 State ကယားပြည်နယ် 9.4
6 7 Kayin MMR003 State ကရင်ပြည်နယ် 9.4
7 8 Magway MMR009 Region မကွေးတိုင်းဒေသကြီး 9.4
8 9 Mandalay MMR010 Region မန္တလေးတိုင်းဒေသကြီး 9.4
9 10 Mon MMR011 State မွန်ပြည်နယ် 9.4
10 11 Nay Pyi Taw MMR018 Union Territory နေပြည်တော် 9.4
geometry
1 MULTIPOLYGON (((95.20798 15...
2 MULTIPOLYGON (((96.17964 19...
3 MULTIPOLYGON (((93.36931 24...
4 MULTIPOLYGON (((97.59674 28...
5 MULTIPOLYGON (((97.1759 19....
6 MULTIPOLYGON (((97.81508 16...
7 MULTIPOLYGON (((94.11699 22...
8 MULTIPOLYGON (((96.14023 23...
9 MULTIPOLYGON (((97.73689 15...
10 MULTIPOLYGON (((96.32013 20...M_District_Sf <- st_read(dsn="data/raw/geospatial", layer = "mmr_polbnda_adm2_250k_mimu") Reading layer `mmr_polbnda_adm2_250k_mimu' from data source
`C:\Users\jiale\Desktop\IS415\IS415-GAA\Take_Home_Exercises\Take_Home_Exercise_1\data\raw\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84M_District_SfSimple feature collection with 80 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 92.1721 ymin: 9.696844 xmax: 101.17 ymax: 28.54554
Geodetic CRS: WGS 84
First 10 features:
OBJECTID ST ST_PCODE DT DT_PCODE DT_MMR PCode_V
1 1 Ayeyarwady MMR017 Hinthada MMR017D002 ဟင်္သာတခရိုင် 9.4
2 2 Ayeyarwady MMR017 Labutta MMR017D004 လပွတ္တာခရိုင် 9.4
3 3 Ayeyarwady MMR017 Maubin MMR017D005 မအူပင်ခရိုင် 9.4
4 4 Ayeyarwady MMR017 Myaungmya MMR017D003 မြောင်းမြခရိုင် 9.4
5 5 Ayeyarwady MMR017 Pathein MMR017D001 ပုသိမ်ခရိုင် 9.4
6 6 Ayeyarwady MMR017 Pyapon MMR017D006 ဖျာပုံခရိုင် 9.4
7 7 Bago (East) MMR007 Bago MMR007D001 ပဲခူးခရိုင် 9.4
8 8 Bago (East) MMR007 Taungoo MMR007D002 တောင်ငူခရိုင် 9.4
9 9 Bago (West) MMR008 Pyay MMR008D001 ပြည်ခရိုင် 9.4
10 10 Bago (West) MMR008 Thayarwady MMR008D002 သာယာဝတီခရိုင် 9.4
geometry
1 MULTIPOLYGON (((95.12637 18...
2 MULTIPOLYGON (((95.04462 15...
3 MULTIPOLYGON (((95.38231 17...
4 MULTIPOLYGON (((94.6942 16....
5 MULTIPOLYGON (((94.27572 15...
6 MULTIPOLYGON (((95.20798 15...
7 MULTIPOLYGON (((95.90674 18...
8 MULTIPOLYGON (((96.17964 19...
9 MULTIPOLYGON (((95.70458 19...
10 MULTIPOLYGON (((95.85173 18...2.2.2.1 Understanding the data set fields.
| Field Name | Description |
| OBJECTID | This is a unique identifier for each feature in the dataset, typically used to identify individual records or polygons in the shapefile. |
| ST | This represents the State or Region in Myanmar. For example, in your dataset, “Ayeyarwady” refers to a state/region. |
| ST_PCODE | This stands for State Postal Code. It is a standardized code that represents each state or region in Myanmar, such as “MMR017” for Ayeyarwady. |
| DT | This stands for District or Township within the respective state/region. For example, “Hinthada” is a district or township within Ayeyarwady. |
| DT_PCODE | This stands for District/Township Postal Code. It is a standardized postal code for each district or township, such as “MMR017D002” for the Hinthada district/township in Ayeyarwady. |
| DT_MMR | This field could be the District/Township name in Myanmar script, written in the local language. It may be an alternative representation of the “DT” field, showing the name of the district/township in Myanmar’s native language. |
| PCODE_V | This could be a Version of the Postal Code or a verification value used internally in the dataset. In this case, the value is “9.4”, possibly indicating a specific version of postal codes or an accuracy measure. |
| geometry | This column represents the spatial data for each district/township. It contains the geometrical shape (MULTIPOLYGON) defining the boundaries of the state or district/township, with coordinates provided in longitude and latitude. |
3.0 Data Pre-processing
To ensure accuracy and usability of the data, several preprocessing steps will be undertaken for the different datasets.
3.1 Myanmar Shapefile
3.1.1 Setting the CRS for the
Since Myanmar uses CRS of 32647 and when we download the map it’s in WGS84, we should change it to 32647 .
st_crs(M_State_Sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_State_Sf <- st_transform(M_State_Sf, crs = 32637)
# Verify that the CRS has been correctly set
print(st_crs(M_State_Sf))Coordinate Reference System:
User input: EPSG:32637
wkt:
PROJCRS["WGS 84 / UTM zone 37N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 37N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",39,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 36°E and 42°E, northern hemisphere between equator and 84°N, onshore and offshore. Djibouti. Egypt. Eritrea. Ethiopia. Georgia. Iraq. Jordan. Kenya. Lebanon. Russian Federation. Saudi Arabia. Somalia. Sudan. Syria. Türkiye (Turkey). Ukraine."],
BBOX[0,36,84,42]],
ID["EPSG",32637]]st_crs(M_District_Sf)Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["latitude",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["longitude",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]]# Set the CRS for m_sf, assuming the appropriate CRS is WGS 84 (EPSG:32647)
M_District_Sf <- st_transform(M_District_Sf, crs = 32637)
# Verify that the CRS has been correctly set
print(st_crs(M_District_Sf))Coordinate Reference System:
User input: EPSG:32637
wkt:
PROJCRS["WGS 84 / UTM zone 37N",
BASEGEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["UTM zone 37N",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",39,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Navigation and medium accuracy spatial referencing."],
AREA["Between 36°E and 42°E, northern hemisphere between equator and 84°N, onshore and offshore. Djibouti. Egypt. Eritrea. Ethiopia. Georgia. Iraq. Jordan. Kenya. Lebanon. Russian Federation. Saudi Arabia. Somalia. Sudan. Syria. Türkiye (Turkey). Ukraine."],
BBOX[0,36,84,42]],
ID["EPSG",32637]]3.1.2 Renaming and removal of column names
colnames(M_State_Sf) <- c("OBJECTID", "state","state_code","type", "state_myr", "mimi_version", "geometry")
M_State_Sf_Cleansed <- M_State_Sf %>% select(state, type, state_myr ,geometry)colnames(M_District_Sf) <- c("OBJECTID", "state", "state_code", "district", "district_code", "district_mmr", "mimi_version", "geometry")
M_District_Sf_Cleansed <- M_District_Sf %>% select(state, district, district_mmr, geometry)3.1.3 Checking for validity of data
When working with spatial data, it’s crucial to ensure that all geometries are valid. Invalid geometries can cause errors in analysis and visualization.
Checking Validity with
st_is_valid():Identifying Invalid Geometries:
Fixing Invalid Geometries with
st_make_valid()
#checking if it's valid
M_State_Sf_Validity <- st_is_valid(M_State_Sf_Cleansed)
M_State_Sf_Invalid <- which(!M_State_Sf_Validity)
if (length(M_State_Sf_Invalid) > 0) {
print("MPZ Invalid!")
print(M_State_Sf_cleansed[M_State_Sf_Invalid, ])
} else {
print("it's valid!")
}[1] "it's valid!"#checking if it's valid
M_District_Sf_Validity <- st_is_valid(M_District_Sf_Cleansed)
M_District_Sf_Invalid <- which(!M_District_Sf_Validity)
if (length(M_District_Sf_Invalid) > 0) {
print("MPZ Invalid!")
print(M_District_Sf_cleansed[M_District_Sf_Invalid, ])
} else {
print("it's valid!")
}[1] "it's valid!"3.1.4 Visualizing the mynamar map
On the top right, you can toggle between the district level and also the state level to understand more about the boundaries of Myanmar.
tmap_mode("plot")
tm_shape(M_State_Sf_Cleansed) + # Base map (Myanmar boundaries)
tm_polygons("state", # Color the base map by state
palette = "Set3", # Use Set3 color palette for states
border.col = "gray", # Border color for the states
alpha = 0.5, # Semi-transparent polygons
title = "State", # Legend title for states
legend.show = TRUE, # Show legend for state colors
) +
tm_layout(main.title = "States of Myanmar", # Main map title
legend.outside = TRUE) # Position the legend outside the map
There is more than 80 district here, so it only showcases 30 :)
tm_shape(M_District_Sf_Cleansed) + # Base map (Myanmar boundaries)
tm_polygons("district", # Color the base map by district
palette = "Set3", # Use Set3 color palette for districts
border.col = "gray", # Border color for the districts
alpha = 0.5, # Semi-transparent polygons
title = "District", # Legend title for districts
legend.show = TRUE, # Show legend for district colors
) +
tm_layout(main.title = "Districts of Myanmar", # Main map title
legend.outside = TRUE) # Position the legend outside the mapWarning: Number of levels of the variable "district" is 80, which is larger
than max.categories (which is 30), so levels are combined. Set
tmap_options(max.categories = 80) in the layer function to show all levels.
3.2 ACLED Data
3.2.1 Changing the Column Names
Since Myanmar’s regional hierarchy follows State, District, and Township levels, we will rename the columns accordingly:
admin1→ Stateadmin2→ Districtadmin3→ Township
This is important because different countries use different administrative hierarchies. For example, in Singapore, the hierarchy is organized by Region and Subzones. Adjusting these names ensures that our dataset aligns with Myanmar’s specific regional structure for accurate analysis.
ACLEDData_Cleanse <- ACLEDData %>%
select(event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1,
actor2, inter2, interaction, admin1, admin2, admin3, location, latitude,
longitude, fatalities) %>%
rename(state = admin1, district = admin2, township = admin3) %>%
mutate(across(where(is.character), ~ replace_na(.x, "NA")), # Replace NA in character columns with "NA"
across(where(is.numeric), ~ replace_na(.x, 0))) # Replace NA in numeric columns with 03.2.2 Adding a “Quarter-Year” Column
To facilitate our temporal analysis, we need to add a “quarter-year” column based on the event_date field. This can be done by adjusting the date format to represent the quarter and year, ensuring that each event is categorized by the specific quarter it occurred in (e.g., Q1-2021, Q2-2022). This will allow for easier analysis of conflict trends over time.
# Convert event_date to Date format (if it's not already a date)
ACLEDData_Cleanse$event_date <- as.Date(ACLEDData_Cleanse$event_date, format="%d-%b-%y") # Adjust the format if needed
# Add a new column that shows the quarter and year
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
mutate(quarter_year = paste0("Q", quarter(event_date), "-", year(event_date)))
head(ACLEDData_Cleanse)# A tibble: 6 × 18
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political viol… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
# ℹ 10 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
# district <chr>, township <chr>, location <chr>, latitude <dbl>,
# longitude <dbl>, fatalities <dbl>, quarter_year <chr>3.2.3 Joining ACLED’s Codebook Description
ACLED’s stores their data for the column “interaction” and “inter1” and “inter2” in codes, using their code book, let’s reorganise their data for simplier view, we can reference the code book here to know what each code represent. Map it out as a csv file read it in and change accordingly.
3.2.3.1 Left joining inter1 and inter’s description.
For more details about each inter code read here.
ACLEDActorInterCode <- read_csv("data/raw/aspatial/ActorTypesInterCode.csv")Rows: 8 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
left_join(ACLEDActorInterCode, by = c("inter1" = "code")) %>%
rename(inter1_description = Description)
# Left join again for inter2
ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
left_join(ACLEDActorInterCode, by = c("inter2" = "code")) %>%
rename(inter2_description = Description)
datatable(head(ACLEDData_Cleanse, 5), options = list(pageLength = 5, autoWidth = TRUE))3.2.3.2 Left joining interaction description.
For more details about each interaction code read here.
ACLEDInteractionCode <- read_csv("data/raw/aspatial/AcledInteractionCodes.csv")Rows: 44 Columns: 2
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Description
dbl (1): code
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.ACLEDData_Cleanse <- ACLEDData_Cleanse %>%
left_join(ACLEDInteractionCode, by = c("interaction" = "code")) %>%
rename(interaction_description = Description)
head(ACLEDData_Cleanse)# A tibble: 6 × 21
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political viol… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political viol… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic deve… Strategic… Milit… 1 NA
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, state <chr>,
# district <chr>, township <chr>, location <chr>, latitude <dbl>,
# longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
# inter1_description <chr>, inter2_description <chr>,
# interaction_description <chr>3.2.3 Making it a SF Object and reverse geolocate state and district
Since ACLED provides longitude and latitude data, I prefer to reverse geolocate the points to match Myanmar’s official state and district boundaries. We are uncertain how ACLED assigns these regions, so to ensure consistency, we remove the original state and district columns from the ACLED data and replace them with the geolocated values.
Steps:
Convert ACLED Data to an SF Object: Using longitude and latitude coordinates, transform the ACLED dataset into a spatial object. Remember taht we have to set CRS 32647 here as well.
Perform a Spatial Join: Match the points from ACLED with the corresponding state and district boundaries from the
m_sfshapefile, selecting only those columns.Remove Original Columns: After the spatial join, drop the original
state,district, andtownshipcolumns from the ACLED dataset.Rename the Joined Columns: Rename the newly joined
state.yanddistrict.ytostateanddistrict, effectively replacing the original columns with the reverse-geolocated values.
# Step 1: Convert ACLEDDataCleanse to an sf object using longitude and latitude
ACLEDData_Cleanse_Sf <- st_as_sf(ACLEDData_Cleanse, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
# Step 2: Transform ACLEDDataCleanse_sf to the same CRS as m_sf (EPSG: 32637)
ACLEDData_Cleanse_Sf <- st_transform(ACLEDData_Cleanse_Sf, crs = 32637)
# Step 3: Perform a spatial join, selecting only the state and district from m_sf
reverse_geolocated_state <- st_join(ACLEDData_Cleanse_Sf, M_State_Sf_Cleansed[, c("state")], join = st_intersects)
# Step 2: Spatial join to add 'district' from M_District_Sf_Cleansed
reverse_geolocated_district <- st_join(reverse_geolocated_state, M_District_Sf_Cleansed[, c("district")], join = st_intersects)
# Step 3: Remove original 'state', 'district', and 'township' columns from ACLEDData_Cleanse (if they exist)
# This step removes the original columns, and then renames the newly joined columns
ACLEDData_Cleanse_Sf <- reverse_geolocated_district %>%
select(-contains("state.x"), -contains("district.x"), -contains("township")) %>% # Remove old state, district, township columns
rename(state = state.y, district = district.y) # Rename newly joined columns
ACLEDData_Cleanse <- st_drop_geometry(ACLEDData_Cleanse_Sf)
# View the updated data
print(ACLEDData_Cleanse_Sf)Simple feature collection with 42608 features and 20 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 6629596 ymin: 2135072 xmax: 8594728 ymax: 5015434
Projected CRS: WGS 84 / UTM zone 37N
# A tibble: 42,608 × 21
event_id_cnty event_date year disorder_type event_type actor1 inter1 actor2
<chr> <date> <dbl> <chr> <chr> <chr> <dbl> <chr>
1 MMR64313 2024-06-30 2024 Political vio… Battles Peopl… 3 Milit…
2 MMR64320 2024-06-30 2024 Political vio… Battles Peopl… 3 Milit…
3 MMR64321 2024-06-30 2024 Political vio… Battles Peopl… 3 Milit…
4 MMR64322 2024-06-30 2024 Strategic dev… Strategic… Milit… 1 NA
5 MMR64323 2024-06-30 2024 Political vio… Battles PKDF … 3 Milit…
6 MMR64324 2024-06-30 2024 Strategic dev… Strategic… Milit… 1 NA
7 MMR64325 2024-06-30 2024 Political vio… Battles Milit… 1 PSLF/…
8 MMR64326 2024-06-30 2024 Political vio… Battles PSLF/… 2 Milit…
9 MMR64328 2024-06-30 2024 Political vio… Battles Milit… 1 PSLF/…
10 MMR64330 2024-06-30 2024 Political vio… Battles Milit… 1 PSLF/…
# ℹ 42,598 more rows
# ℹ 13 more variables: inter2 <dbl>, interaction <dbl>, location <chr>,
# latitude <dbl>, longitude <dbl>, fatalities <dbl>, quarter_year <chr>,
# inter1_description <chr>, inter2_description <chr>,
# interaction_description <chr>, geometry <POINT [m]>, state <chr>,
# district <chr>3.2.5 Visualizing it by Event Type
# Set tmap mode to plot for static maps
tmap_mode("plot")tmap mode set to plottingbbox <- st_bbox(M_State_Sf_Cleansed) # Assuming you are using sf object
# Create the tmap object with the base map and event markers
tm_shape(M_State_Sf_Cleansed) +
tm_polygons("state", alpha = 0.3, border.col = "gray") + # Base map with higher transparency
# Add event markers (bubbles) with size based on fatalities
tm_shape(ACLEDData_Cleanse_Sf) +
tm_bubbles(size = "fatalities", # Marker size based on fatalities
col = "event_type", # Color markers by event type
palette = "Set1", # Use Set1 color palette for event types
border.col = "black", # Border color for bubbles
border.alpha = 0.5, # Semi-transparent border
title.size = "Number of Fatalities", # Title for bubble size legend
title.col = "Event Types") + # Title for event type legend
# Layout settings for the map, including title and legend positioning
tm_layout(main.title = "Myanmar's State Conflicts by Fatalities", # Main map title
main.title.size = 1.5, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.2, # Size of the legend title
legend.text.size = 1, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05)) 
4.0 Exploratory Data analysis
With the data cleansed, let’s conduct a high-level analysis of our dataset to determine some key statistics.
The central question we’re exploring is:
“Which district in Myanmar has the highest number of conflicts that may be concerning to civilians?”
In section 3.2.5, we recognized that Myanmar has many states and numerous conflicts during this period. Therefore, we will aggregate the data at the state level to prioritize states based on the following key figures:
Civilian Involvement: We focus on events where either Actor 1 or Actor 2 is a civilian, as these are likely to raise humanitarian concerns.
Total Conflicts: This measures the overall number of conflict events in each state, helping us identify which states experienced the most conflict.
Conflict Density: This calculates the number of conflicts per square kilometer in each state, allowing us to understand the intensity of conflicts relative to the size of the state.
Total Fatalities: This sums up the total number of fatalities within each state, highlighting areas with the highest death toll, which are likely experiencing the most severe impacts.
Fatalities Density: This measures the fatalities per square kilometer, giving a sense of how deadly the conflicts are in each state in relation to its area.
This analysis will allow us to focus on the states most affected by conflicts that pose the greatest risk to civilians.
4.1 Aggregation of Data for exploratory purposes
Step 1: Filter Data for Civilians
In this first step, the data (ACLEDData_Cleanse) is filtered so that only rows where either “Civilians” are involved in the conflict (inter1_description or inter2_description) are kept. Additionally, rows where the state information is missing (NA) are removed because these conflicts can’t be assigned to a specific geographical area.
filtered_data <- ACLEDData_Cleanse %>%
filter((inter1_description == "Civilians" | inter2_description == "Civilians") & !is.na(state))Step 2: Calculate State Area in Square Kilometers
Next, the code calculates the area for each state from the M_State_Sf_Cleansed spatial dataset. The st_area function retrieves the area in square meters, and dividing by 1 million (1e6) converts it into square kilometers. This information will be used later to calculate conflict and fatality density.
state_area_km2 <- st_area(M_State_Sf_Cleansed) / 1e6
state_area_df <- data.frame(state = M_State_Sf_Cleansed$state, area_km2 = as.numeric(state_area_km2))Step 3: Aggregate Total Conflicts and Fatalities by State
The next step is to group the filtered data by state and calculate two key metrics: the total number of conflicts and the total number of fatalities for each state. This is done using group_by and summarise.
agg_data_by_state <- filtered_data %>%
group_by(state) %>%
summarise(
total_conflicts = n(),
total_fatalities = sum(fatalities, na.rm = TRUE),
.groups = "drop"
)Step 4: Convert to Non-Spatial Data Frame
Here, the spatial data is converted into a regular data frame by dropping the geometry information. This allows easier manipulation of the non-spatial data for further processing.
filtered_data_df <- st_drop_geometry(filtered_data)Step 5: Prepare Data Frame for Yearly Data
To store yearly conflict and fatality data for each state, an empty data frame is initialized by selecting distinct states from the filtered data
final_result <- filtered_data_df %>% select(state) %>% distinct()Step 6: Loop Through Each Year and Calculate Yearly Conflicts/Fatalities
In this step, the code iterates over each unique year in the dataset. For each year, it calculates the total conflicts and fatalities for each state. The results are merged back into the final_result data frame, with columns named according to the year (e.g., conflicts_2021, fatalities_2021).
unique_years <- unique(filtered_data_df$year)
for (yr in unique_years) {
yearly_data <- filtered_data_df %>%
filter(year == yr) %>%
group_by(state) %>%
summarise(
conflicts = n(),
fatalities = sum(fatalities, na.rm = TRUE)
) %>%
rename_at(vars(conflicts, fatalities), ~ paste0(., "_", yr)) # Rename with year
final_result <- left_join(final_result, yearly_data, by = "state")
}Step 7: Replace NAs with 0
Any missing data for a state-year combination (e.g., a state with no conflicts in a particular year) is replaced with 0 to avoid leaving gaps in the analysis.
final_result <- final_result %>%
mutate(across(everything(), ~ ifelse(is.na(.), 0, .)))Step 8: Calculate Conflict and Fatality Density
Now, the code calculates the density of conflicts and fatalities per square kilometer for each state. The results from Step 3 (agg_data_by_state) are joined with the state area data from Step 2. The density is simply the number of conflicts or fatalities divided by the area in square kilometers.
agg_data_with_density <- agg_data_by_state %>%
left_join(state_area_df, by = "state") %>%
mutate(
conflict_density = total_conflicts / area_km2,
fatality_density = total_fatalities / area_km2
)Step 9: Merge Yearly Data and Density with Spatial Data
Finally, the spatial data (M_State_Sf_Cleansed) is merged with the yearly conflict/fatality data (final_result) and the density data (agg_data_with_density). Additional calculations are made to find the fatality-per-conflict ratio which are then added to the final spatial dataset.
finalized_map <- M_State_Sf_Cleansed %>%
left_join(final_result, by = "state") %>%
left_join(agg_data_with_density, by = "state")
finalized_map <- finalized_map %>%
mutate(conflict_fatality_ratio = ifelse(total_fatalities == 0, NA, total_fatalities / total_conflicts))
finalized_map <- finalized_map %>%
mutate(
fatality_per_conflict = ifelse(total_conflicts == 0, NA, total_fatalities / total_conflicts),
)
finalized_map <- finalized_map %>%
mutate(
rank_total_conflicts = rank(-total_conflicts, ties.method = "min"), # Rank total_conflicts (largest to smallest)
rank_total_fatalities = rank(-total_fatalities, ties.method = "min"), # Rank total_fatalities (largest to smallest)
rank_conflict_density = rank(-conflict_density, ties.method = "min"), # Rank conflict_density (largest to smallest)
rank_fatality_density = rank(-fatality_density, ties.method = "min"), # Rank fatality_density (largest to smallest)
rank_fatality_per_conflict = rank(-fatality_per_conflict, ties.method = "min") # Rank fatality_per_conflict (largest to smallest)
)Final Output: Display the Data
The finalized_map contains all the processed information, including conflict/fatality counts, densities, and state geometry, which can now be visualized or further analyzed.
head(finalized_map)Simple feature collection with 6 features and 23 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 6589238 ymin: 2946007 xmax: 8116377 ymax: 5139186
Projected CRS: WGS 84 / UTM zone 37N
state type state_myr conflicts_2024 fatalities_2024
1 Ayeyarwady Region ဧရာဝတီတိုင်းဒေသကြီး 64 9
2 Bago Region ပဲခူးတိုင်းဒေသကြီး 194 128
3 Chin State ချင်းပြည်နယ် 72 52
4 Kachin State ကချင်ပြည်နယ် 160 95
5 Kayah State ကယားပြည်နယ် 21 38
6 Kayin State ကရင်ပြည်နယ် 52 59
conflicts_2023 fatalities_2023 conflicts_2022 fatalities_2022 conflicts_2021
1 83 7 177 29 233
2 339 220 172 61 191
3 101 76 204 80 160
4 316 133 226 140 221
5 113 74 149 87 114
6 165 83 165 86 116
fatalities_2021 total_conflicts total_fatalities area_km2 conflict_density
1 34 557 79 91740.37 0.006071482
2 87 896 496 106296.56 0.008429247
3 32 537 240 83773.78 0.006410120
4 71 923 439 214569.60 0.004301635
5 108 397 307 33443.39 0.011870806
6 57 498 285 92162.25 0.005403514
fatality_density geometry conflict_fatality_ratio
1 0.0008611258 MULTIPOLYGON (((7515288 297... 0.1418312
2 0.0046661904 MULTIPOLYGON (((7376294 368... 0.5535714
3 0.0028648583 MULTIPOLYGON (((6594137 416... 0.4469274
4 0.0020459562 MULTIPOLYGON (((6707242 513... 0.4756230
5 0.0091796915 MULTIPOLYGON (((7484860 384... 0.7732997
6 0.0030923725 MULTIPOLYGON (((7888506 331... 0.5722892
fatality_per_conflict rank_total_conflicts rank_total_fatalities
1 0.1418312 11 14
2 0.5535714 9 6
3 0.4469274 12 13
4 0.4756230 8 7
5 0.7732997 14 10
6 0.5722892 13 12
rank_conflict_density rank_fatality_density rank_fatality_per_conflict
1 11 15 15
2 8 8 4
3 10 10 9
4 14 12 6
5 7 4 1
6 12 9 34.2 Visualizing the aggregated data.
For each column that we created earlier, we can now showcase it in a choropleth map to highlight the states with the highest values in each category. These maps will help visually identify which state has the highest value for each column. You can navigate through the tabset to explore the different metrics and easily view the data on a map!
# Create the tmap object for total civilian conflicts
tm_total_conflicts <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
# Layer for total civilian conflicts
tm_shape(finalized_map) +
tm_polygons("total_conflicts",
title = "Total Civilian Conflicts",
palette = "Reds",
border.col = "black",
alpha = 0.8, # Transparency for polygons
border.alpha = 0.5, # Semi-transparent border
id = "state") + # Set state as identifier for the polygons
# Layout settings matching your style
tm_layout(
main.title = "Myanmar's State by Total Civilian Related Conflicts", # Main map title
main.title.size = 1.5, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.2, # Size of the legend title
legend.text.size = 1, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05) # Set inner margins for better spacing
)
tm_total_conflicts
# Create the tmap object for total civilian conflicts
tm_total_fatalities <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
# Layer for total civilian conflicts
tm_shape(finalized_map) +
tm_polygons("total_fatalities",
title = "Total Fatalities by Civilian Related Conflicts",
palette = "Blues",
border.col = "black",
alpha = 0.8, # Transparency for polygons
border.alpha = 0.5, # Semi-transparent border
id = "state") + # Set state as identifier for the polygons
# Layout settings matching your style
tm_layout(
main.title = "Myanmar's State by Total Fatatlies For Civilian Related Conflicts", # Main map title
main.title.size = 1.0, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.2, # Size of the legend title
legend.text.size = 1, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05) # Set inner margins for better spacing
)
tm_total_fatalities
tm_conflict_density <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
# Layer for total civilian conflicts
tm_shape(finalized_map) +
tm_polygons("conflict_density",
title = "Myanmar State by Civilian Related Conflicts Density",
palette = "Blues",
border.col = "black",
alpha = 0.8, # Transparency for polygons
border.alpha = 0.5, # Semi-transparent border
id = "state") + # Set state as identifier for the polygons
# Layout settings matching your style
tm_layout(
main.title = "Myanmar State by Civilian Related Conflicts Density", # Main map title
main.title.size = 1.0, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.2, # Size of the legend title
legend.text.size = 1, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05) # Set inner margins for better spacing
)
tm_conflict_density
tm_fatalities_conflict <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
# Layer for total civilian conflicts
tm_shape(finalized_map) +
tm_polygons("fatality_density",
title = "Myanmar State by Fatalies For Civilians Related Conflicts Density Per Km^2",
palette = "Blues",
border.col = "black",
alpha = 0.8, # Transparency for polygons
border.alpha = 0.5, # Semi-transparent border
id = "state") + # Set state as identifier for the polygons
# Layout settings matching your style
tm_layout(
main.title = "Myanmar State by Fatalies For Civilians Related Conflicts Density Per Km^2", # Main map title
main.title.size = 1.0, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.0, # Size of the legend title
legend.text.size = 0.8, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05) # Set inner margins for better spacing
)
tm_fatalities_conflict
tm_fatalities_density <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
# Layer for total civilian conflicts
tm_shape(finalized_map) +
tm_polygons("fatality_per_conflict",
title = "Myanmar State by Fatalies for Civilian Related Conflicts",
palette = "Blues",
border.col = "black",
alpha = 0.8, # Transparency for polygons
border.alpha = 0.5, # Semi-transparent border
id = "state") + # Set state as identifier for the polygons
# Layout settings matching your style
tm_layout(
main.title = "Myanmar State by Fatalies for Civilian Related Conflicts", # Main map title
main.title.size = 1.0, # Title font size
legend.outside = TRUE, # Position legend outside the map
legend.outside.size = 0.5, # Adjust size of outside legend
legend.position = c("left", "top"), # Position for the event type legend
legend.title.size = 1.0, # Size of the legend title
legend.text.size = 0.8, # Size of the legend text
legend.bg.color = "white", # Background color for the legend
legend.bg.alpha = 0.5, # Transparency for the legend background
inner.margins = c(0.05, 0.05, 0.05, 0.05) # Set inner margins for better spacing
)
tm_fatalities_density
4.3 Analysis on the aggregated data to find top 3 state to analyse
4.3.1 How to choose which States?
With the maps above, everyone can choose what they would like to analyze, whether it’s total civilian conflicts or fatalities by state. However, I would like to focus on the central question:
Which state has the highest amount of conflict and fatality density?
By using density, I’ve aimed to assess fatalities in relation to conflicts and ensure a fair comparison by adjusting for density.
Note
Keep in mind, this is not the most precise analysis, as other factors—such as ethnicity, population demographics, and gender—could also be considered. However, for the purpose of this study, I have chosen to focus on conflict and fatality density.
Conflict Density: This is the number of conflicts per square kilometer in a state.
Conflict Density=Total ConflictsArea of State (km²) = Conflict Density=Area of State (km²)Total Conflicts
This tells you how many conflicts occur per unit of area (1 km²) in each state.
Fatality Density: This is the number of fatalities per square kilometer.
Fatality Density=Total FatalitiesArea of State (km²) = Fatality Density=Area of State (km²)Total Fatalities
4.3.2 Viewing the top 3 states across
Using the map below, you can interactively explore the Conflict and Fatality Density. We observe that the top three states are as follows, and here’s how they compare across the data set.
You can use the layer toggle to switch between fatalities and conflicts, or click on the map layers to view detailed, aggregated data for each state.
# Step 1: Set tmap mode to view for interactive maps
tmap_mode("view")tmap mode set to interactive viewing# Step 2: Define popup variables to show all relevant columns when clicking on a state
popup_variables <- c(
"Total Civilian Conflicts" = "total_conflicts",
"Fatalities" = "total_fatalities",
"Civilian Conflict Density" = "conflict_density",
"Fatality Density" = "fatality_density",
"Fatality Per Civilian Conflict" = "fatality_per_conflict",
"Ranked For Total Conflict By State (Out of 15)" = "rank_total_conflicts",
"Ranked For Total Fatalities By State (Out of 15)" = "rank_total_fatalities",
"Ranked For Conflict Density By State (Out of 15)" = "rank_conflict_density",
"Ranked For Fatalities Density By State (Out of 15)" = "rank_fatality_density",
"Ranked For Fatalies/Conflict By State (Out of 15)" = "rank_fatality_per_conflict"
)
# Step 3: Create the map with 6 layers for toggling
tm <- tm_basemap(server = c("Esri.WorldGrayCanvas", "OpenStreetMap", "Esri.WorldTopoMap")) +
tm_shape(finalized_map) +
tm_polygons("conflict_density",
title = "Civilian Conflict Density",
palette = "Oranges",
border.col = "black",
popup.vars = popup_variables,
id = "state",
group = "Civilian Conflict Density") +
tm_shape(finalized_map) +
tm_polygons("fatality_density",
title = "Fatality Density",
palette = "Purples",
border.col = "black",
popup.vars = popup_variables,
id = "state",
group = "Fatality Density") +
tm_layout(
legend.outside = TRUE,
legend.outside.size = 0.5, # Adjust the size of the legend
legend.position = c("left", "top") # Position the legend
)
# Step 5: Display the interactive map
tmlegend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.tmap_mode("plot")tmap mode set to plottingfinalized_data_df <- as.data.frame(st_drop_geometry(finalized_map))
# Sort the data by conflict_density and fatality_density in descending order
finalized_data_df <- finalized_data_df %>%
arrange(desc(conflict_density), desc(fatality_density))
# Format the densities to show "per km²"
finalized_data_df$conflict_density_label <- paste0(round(finalized_data_df$conflict_density, 2), " per km²")
finalized_data_df$fatality_density_label <- paste0(round(finalized_data_df$fatality_density, 2), " per km²")
# Create bar chart for Conflict Density (sorted by conflict_density)
conflict_density_plot <- ggplot(finalized_data_df, aes(x = reorder(state, -conflict_density), y = conflict_density)) +
geom_bar(stat = "identity", fill = "orange") +
geom_text(aes(label = conflict_density_label), vjust = 1.8, size = 1.5, color = "black", position = position_stack(vjust = 0.5)) + # Place text at bottom of bar
ggtitle("Top States by Conflict Density (per km²)") +
xlab("State") +
ylab("Conflict Density (per km²)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Create bar chart for Fatality Density (sorted by fatality_density)
fatality_density_plot <- ggplot(finalized_data_df, aes(x = reorder(state, -fatality_density), y = fatality_density)) +
geom_bar(stat = "identity", fill = "purple") +
geom_text(aes(label = fatality_density_label), vjust = 1.8, size = 1.5, color = "black", position = position_stack(vjust = 0.5)) + # Place text at bottom of bar
ggtitle("Top States by Fatality Density (per km²)") +
xlab("State") +
ylab("Fatality Density (per km²)") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Display both charts side by side
grid.arrange(conflict_density_plot, fatality_density_plot, ncol = 2)
4.3.3 Top 3 Chosen state for further analysis
4.3.3.1 Yangon
4.3.3.2 Mandalay
4.3.3.3 Sagaling.
4.4 Final Round of Cleansing ACLED Data
Now that we have finalized the data and narrowed down the states that we are going to analyse, we will filter out for the states data.
Yangon_ACLED_Data <- filtered_data %>%
filter(state %in% c("Yangon"))
class(Yangon_ACLED_Data)[1] "tbl_df" "tbl" "data.frame"Yangon_ACLED_Data_Sf <- st_as_sf(Yangon_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
Yangon_ACLED_Data_Sf <- st_transform(Yangon_ACLED_Data_Sf, crs = 32637)
class(Yangon_ACLED_Data_Sf)[1] "sf" "tbl_df" "tbl" "data.frame"Yangon_ACLED_Validity <- st_is_valid(Yangon_ACLED_Data_Sf)
Yangon_invalid <- which(!Yangon_ACLED_Validity)
if (length(Yangon_invalid) > 0) {
print("Yangon is Invalid!")
print(Yangon_ACLED_Data_Sf[Yangon_invalid, ])
} else {
print("Yangon_ACLED_Data_Sf is valid!")
}[1] "Yangon_ACLED_Data_Sf is valid!"Yangon_ACLED_SFDF <- as_Spatial(Yangon_ACLED_Data_Sf)
Yangon_ACLED_SFDFclass : SpatialPointsDataFrame
features : 1489
extent : 7470244, 7680669, 3162567, 3345642 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs
variables : 20
names : event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1, actor2, inter2, interaction, location, latitude, longitude, fatalities, quarter_year, ...
min values : MMR10949, 18633, 2021, Political violence, Explosions/Remote violence, 5 Brothers Younger, 1, Civilians (Australia), 0, 17, Ah Hpyauk, 16.4389, 95.6934, 0, Q1-2021, ...
max values : MMR64089, 19895, 2024, Strategic developments, Violence against civilians, YUG: Yangon Urban Guerrillas, 7, NA, 7, 70, Zee Hpyu Kone, 17.5612, 96.7448, 13, Q4-2023, ... Mandalay_ACLED_Data <- filtered_data %>%
filter(state %in% c("Mandalay"))
class(Mandalay_ACLED_Data)[1] "tbl_df" "tbl" "data.frame"Mandalay_ACLED_Data_Sf <- st_as_sf(Mandalay_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
class(Mandalay_ACLED_Data_Sf)[1] "sf" "tbl_df" "tbl" "data.frame"Mandalay_ACLED_Validity <- st_is_valid(Mandalay_ACLED_Data_Sf)
Mandalay_invalid <- which(!Mandalay_ACLED_Validity)
if (length(Mandalay_invalid) > 0) {
print("Mandalay is Invalid!")
print(Mandalay_ACLED_Data_Sf[Mandalay_invalid, ])
} else {
print("Mandalay_ACLED_Data_Sf is valid!")
}[1] "Mandalay_ACLED_Data_Sf is valid!"Mandalay_ACLED_SFDF <- as_Spatial(Mandalay_ACLED_Data_Sf)
Mandalay_ACLED_SFDFclass : SpatialPointsDataFrame
features : 2021
extent : 94.849, 96.6148, 20.4319, 23.667 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 20
names : event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1, actor2, inter2, interaction, location, latitude, longitude, fatalities, quarter_year, ...
min values : MMR10960, 18640, 2021, Political violence, Explosions/Remote violence, 21KPG: 21 Guerrilla Force - Kyaukpadaung, 1, Civilians (Myanmar), 0, 17, 4 Maing Kan Thar, 20.4319, 94.849, 0, Q1-2021, ...
max values : MMR64318, 19904, 2024, Strategic developments, Violence against civilians, Zero Guerrilla Force - Myingyan, 7, NA, 7, 70, Zin Chaung, 23.667, 96.6148, 10, Q4-2023, ... Sagaing_ACLED_Data <- filtered_data %>%
filter(state %in% c("Sagaing"))
class(Sagaing_ACLED_Data)[1] "tbl_df" "tbl" "data.frame"Sagaing_ACLED_Data_Sf <- st_as_sf(Sagaing_ACLED_Data, coords = c("longitude", "latitude"), crs = 4326, remove = FALSE)
class(Sagaing_ACLED_Data_Sf)[1] "sf" "tbl_df" "tbl" "data.frame"Sagaing_ACLED_Validity <- st_is_valid(Sagaing_ACLED_Data_Sf)
Sagaing_invalid <- which(!Sagaing_ACLED_Validity)
if (length(Sagaing_invalid) > 0) {
print("Sagaing is Invalid!")
print(Sagaing_ACLED_Data_Sf[Sagaing_invalid, ])
} else {
print("Sagaing_ACLED_Data_Sf is valid!")
}[1] "Sagaing_ACLED_Data_Sf is valid!"Sagaing_ACLED_SFDF <- as_Spatial(Sagaing_ACLED_Data_Sf)
Sagaing_ACLED_SFDFclass : SpatialPointsDataFrame
features : 5346
extent : 93.9575, 96.6034, 21.605, 26.3006 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 20
names : event_id_cnty, event_date, year, disorder_type, event_type, actor1, inter1, actor2, inter2, interaction, location, latitude, longitude, fatalities, quarter_year, ...
min values : MMR10952, 18640, 2021, Political violence, Explosions/Remote violence, ABSDF: All Burma Students' Democratic Front, 1, Civilians (China), 0, 17, 55 Maing, 21.605, 93.9575, 0, Q1-2021, ...
max values : MMR65891, 19898, 2024, Strategic developments, Violence against civilians, Zero Guerrilla Force - Myingyan, 7, NA, 7, 70, Zin Ka Le (Zin Ka Lin), 26.3006, 96.6034, 175, Q4-2023, ... Viola now we have the final state that that en
5.0 Spatial Data Wrangling
Now that we have chosen our 3 states to analyse we have to prepare the data to analyse them within
5.1 Setting up the Owin Window for States.
In this part of the code, we are setting up a spatial window for Myanmar, which defines the boundaries within which the spatial point pattern analysis will take place. This is done using the owin function from the spatstat package, which converts the polygonal boundary of Myanmar into a format that can be used for further spatial analysis.
Yangon_Sf <- M_State_Sf %>%
filter(state %in% c("Yangon"))
# Step 1: Extract the individual polygons from the multipolygon
yangon_polygons <- st_cast(Yangon_Sf$geometry, "POLYGON")
yangon_polygons_filtered <- yangon_polygons[-c(1, 2)]
yangon_multipolygon_filtered <- st_combine(yangon_polygons_filtered)
# Add the filtered multipolygon back to the Yangon_Sf object
Yangon_Sf$geometry <- yangon_multipolygon_filtered
Yangon <- as_Spatial(Yangon_Sf)
Yangon_Owin <- as.owin(Yangon_Sf)
plot(Yangon_Owin)
summary(Yangon_Owin)Window: polygonal boundary
6 separate polygons (1 hole)
vertices area relative.area
polygon 1 3832 2.79639e+10 1.00e+00
polygon 2 (hole) 3 -5.16800e+02 -1.85e-08
polygon 3 14 6.05403e+05 2.16e-05
polygon 4 11 4.56071e+05 1.63e-05
polygon 5 20 5.15327e+06 1.84e-04
polygon 6 37 7.35365e+06 2.63e-04
enclosing rectangle: [7458583, 7698916] x [3150830, 3397886] units
(240300 x 247100 units)
Window area = 27977500000 square units
Fraction of frame area: 0.471Mandalay_Sf <- M_State_Sf %>%
filter(state %in% c("Mandalay"))
Mandalay <- as_Spatial(Mandalay_Sf)
Mandalay_Owin <- as.owin(Mandalay_Sf)
plot(Mandalay_Owin)
summary(Mandalay_Owin)Window: polygonal boundary
single connected closed polygon with 5914 vertices
enclosing rectangle: [6990555, 7369406] x [3760343, 4336524] units
(378900 x 576200 units)
Window area = 79199300000 square units
Fraction of frame area: 0.363Sagaing_Sf <- M_State_Sf %>%
filter(state %in% c("Sagaing"))
Sagaing <- as_Spatial(Sagaing_Sf)
Sagaing_Owin <- as.owin(Sagaing_Sf)
plot(Sagaing_Owin)
summary(Sagaing_Owin)Window: polygonal boundary
single connected closed polygon with 5882 vertices
enclosing rectangle: [6600929, 7143229] x [3910303, 4901317] units
(542300 x 991000 units)
Window area = 2.21637e+11 square units
Fraction of frame area: 0.4125.2 Setting up the ACLED Spatial Class
5.2.1 Converting Spatial Class Into Generic PPP Format
# Yangon_ACLED_Data
# Yangon_ACLED_Data_Sf
# Yangon_ACLED_SFDF
Yangon_ACLED_coords <- st_coordinates(Yangon_ACLED_Data_Sf)
# Define the window using the bounding box
Yangon_ACLED_bbox <- st_bbox(Yangon_ACLED_Data_Sf)
Yangon_ACLED_window <- owin(xrange = Yangon_ACLED_bbox[c("xmin", "xmax")], yrange = Yangon_ACLED_bbox[c("ymin", "ymax")])
# Create the ppp object
Yangon_ACLED_ppp <- ppp(x = Yangon_ACLED_coords[, 1], y = Yangon_ACLED_coords[, 2], window = Yangon_ACLED_window)Warning: data contain duplicated points# Check the ppp object
summary(Yangon_ACLED_ppp)Planar point pattern: 1489 points
Average intensity 3.865152e-08 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 9 decimal places
Window: rectangle = [7470244, 7680669] x [3162567, 3345642] units
(210400 x 183100 units)
Window area = 38523700000 square unitsplot(Yangon_ACLED_ppp)
#Mandalay_ACLED_Data
# Mandalay_ACLED_Data_Sf
# Mandalay_ACLED_SFDF
Mandalay_ACLED_coords <- st_coordinates(Mandalay_ACLED_Data_Sf)
# Define the window using the bounding box
Mandalay_ACLED_bbox <- st_bbox(Mandalay_ACLED_Data_Sf)
Mandalay_ACLED_window <- owin(xrange = Mandalay_ACLED_bbox[c("xmin", "xmax")], yrange = Mandalay_ACLED_bbox[c("ymin", "ymax")])
# Create the ppp object
Mandalay_ACLED_ppp <- ppp(x = Mandalay_ACLED_coords[, 1], y = Mandalay_ACLED_coords[, 2], window = Mandalay_ACLED_window)Warning: data contain duplicated points# Check the ppp object
summary(Mandalay_ACLED_ppp)Planar point pattern: 2021 points
Average intensity 353.7831 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 4 decimal places
Window: rectangle = [94.849, 96.6148] x [20.4319, 23.667] units
(1.766 x 3.235 units)
Window area = 5.71254 square unitsplot(Mandalay_ACLED_ppp)
Sagaing_ACLED_coords <- st_coordinates(Sagaing_ACLED_Data_Sf)
# Define the window using the bounding box
Sagaing_ACLED_bbox <- st_bbox(Sagaing_ACLED_Data_Sf)
Sagaing_ACLED_window <- owin(xrange = Sagaing_ACLED_bbox[c("xmin", "xmax")], yrange = Sagaing_ACLED_bbox[c("ymin", "ymax")])
# Create the ppp object
Sagaing_ACLED_ppp <- ppp(x = Sagaing_ACLED_coords[, 1], y = Sagaing_ACLED_coords[, 2], window = Sagaing_ACLED_window)Warning: data contain duplicated points# Check the ppp object
summary(Sagaing_ACLED_ppp)Planar point pattern: 5346 points
Average intensity 430.2932 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 4 decimal places
Window: rectangle = [93.9575, 96.6034] x [21.605, 26.3006] units
(2.646 x 4.696 units)
Window area = 12.4241 square unitsplot(Sagaing_ACLED_ppp)
5.2.2 Dealing with duplicate conflicts
5.2.2.1 Check For Duplicate
any(duplicated(Yangon_ACLED_ppp))[1] TRUEsum(multiplicity(Yangon_ACLED_ppp) > 1)[1] 1432any(duplicated(Yangon_ACLED_ppp))[1] TRUEsum(multiplicity(Yangon_ACLED_ppp) > 1)[1] 14325.2.2.2 Unique Marking them
Yangon_marks <- rep(1, npoints(Yangon_ACLED_ppp))
Yangon_marks[duplicated(Yangon_ACLED_ppp)] <- 2
# Create a new ppp object with marks attached to each point
marked_Yangon_ACLED_ppp <- ppp(x = Yangon_ACLED_coords[, 1], y = Yangon_ACLED_coords[, 2], window = Yangon_ACLED_window, marks = Yangon_marks)Warning: data contain duplicated points# Check the ppp object with unique marks
summary(marked_Yangon_ACLED_ppp)Marked planar point pattern: 1489 points
Average intensity 3.865152e-08 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 9 decimal places
marks are numeric, of type 'double'
Summary:
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 2.000 2.000 1.916 2.000 2.000
Window: rectangle = [7470244, 7680669] x [3162567, 3345642] units
(210400 x 183100 units)
Window area = 38523700000 square unitsplot(marked_Yangon_ACLED_ppp, main = "Yangon Conflicts with Unique Marks", cols = c("black", "red"))
Note
Why use mark and not jitter or delete for these duplicates?
When dealing with conflict data, such as ACLED, uniquely marking events is essential to preserve both the spatial and temporal accuracy of key occurrences. Events often take place at the same location but during different time periods, and jittering the data would distort their true geographic positions, while deleting them would erase important historical information. By marking events as unique, we ensure that the dataset retains its integrity, allowing for the analysis of recurring conflicts without compromising precision. This approach ensures that both spatial and temporal patterns of conflict hotspots are maintained, which is crucial for understanding the dynamics of the events over time.
5.3.1 Combining Point Events and State’s Owin Object
marked_Yangon_ACLED_ppp[Yangon_Owin]Marked planar point pattern: 1489 points
marks are numeric, of storage type 'double'
window: polygonal boundary
enclosing rectangle: [7458583, 7698916] x [3150830, 3397886] unitsplot(marked_Yangon_ACLED_ppp[Yangon_Owin])
5.4 Combining them both
# Step 1: Define the bounding box and window for Yangon
Yangon_ACLED_bbox <- st_bbox(Yangon_Sf)
Yangon_ACLED_window <- owin(
xrange = c(Yangon_ACLED_bbox["xmin"], Yangon_ACLED_bbox["xmax"]),
yrange = c(Yangon_ACLED_bbox["ymin"], Yangon_ACLED_bbox["ymax"])
)
# Extract unique quarters
quarters <- unique(Yangon_ACLED_Data_Sf$quarter_year)
# Initialize an empty list to store ppp objects
plot_list <- list()
# Step 2: Loop over each quarter, jitter duplicates, and store ppp objects
for (quarter in quarters) {
# Filter the data for the current quarter
quarter_data <- Yangon_ACLED_Data_Sf %>%
filter(quarter_year == quarter)
# Extract coordinates for this quarter
quarter_coords <- st_coordinates(quarter_data)
# Create a ppp object for this quarter
quarter_ppp <- ppp(
x = quarter_coords[, 1],
y = quarter_coords[, 2],
window = Yangon_ACLED_window
)
# Step 3: Jitter duplicates
is_duplicate <- duplicated(quarter_ppp)
jittered_x <- quarter_coords[, 1] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
jittered_y <- quarter_coords[, 2] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
# Create a jittered ppp object with the new coordinates
jittered_conflict_ppp <- ppp(
x = jittered_x,
y = jittered_y,
window = Yangon_ACLED_window
)
# Step 4: Subset the jittered points within Yangon window (Yangon owins)
quarter_conflict_ppp <- jittered_conflict_ppp[Yangon_Owin]
# Step 5: Store the ppp object for later plotting
plot_list[[quarter]] <- quarter_conflict_ppp
}Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points# Step 6: Plot all quarter data using spatstat's plot function
par(mfrow = c(4, 4)) # Adjust rows and columns for multiple plots
for (quarter in names(plot_list)) {
plot(plot_list[[quarter]], main = paste("Yangon Conflict Events -", quarter))
}
Yangon_ACLED_bbox <- st_bbox(Yangon_Sf)
Yangon_ACLED_window <- owin(
xrange = c(Yangon_ACLED_bbox["xmin"], Yangon_ACLED_bbox["xmax"]),
yrange = c(Yangon_ACLED_bbox["ymin"], Yangon_ACLED_bbox["ymax"])
)
# Extract unique quarters
quarters <- unique(Yangon_ACLED_Data_Sf$quarter_year)
# Initialize an empty dictionary-like list to store ppp objects by quarter
yangon_quarters_conflict <- list()
# Step 2: Loop over each quarter, jitter duplicates, and store the ppp object in the list
for (quarter in quarters) {
# Filter the data for the current quarter
quarter_data <- Yangon_ACLED_Data_Sf %>%
filter(quarter_year == quarter)
# Extract coordinates for this quarter
quarter_coords <- st_coordinates(quarter_data)
# Create a ppp object for this quarter
quarter_ppp <- ppp(
x = quarter_coords[, 1],
y = quarter_coords[, 2],
window = Yangon_ACLED_window
)
# Step 3: Jitter duplicates
is_duplicate <- duplicated(quarter_ppp)
jittered_x <- quarter_coords[, 1] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
jittered_y <- quarter_coords[, 2] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
# Create a new ppp object with jittered points
jittered_ppp <- ppp(
x = jittered_x,
y = jittered_y,
window = Yangon_ACLED_window
)
# Step 4: Store the jittered ppp object in the dictionary (list) with the quarter as the key
yangon_quarters_conflict[[quarter]] <- jittered_ppp
}Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated pointskde_plot_list <- list()
# Step 2: Loop over each quarter, jitter duplicates, and store KDE objects
for (quarter in quarters) {
# Filter the data for the current quarter
quarter_data <- Yangon_ACLED_Data_Sf %>%
filter(quarter_year == quarter)
# Extract coordinates for this quarter
quarter_coords <- st_coordinates(quarter_data)
# Create a ppp object for this quarter
quarter_ppp <- ppp(
x = quarter_coords[, 1],
y = quarter_coords[, 2],
window = Yangon_ACLED_window
)
# Step 3: Jitter duplicates
is_duplicate <- duplicated(quarter_ppp)
jittered_x <- quarter_coords[, 1] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
jittered_y <- quarter_coords[, 2] + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
# Create a jittered ppp object with the new coordinates
jittered_conflict_ppp <- ppp(
x = jittered_x,
y = jittered_y,
window = Yangon_ACLED_window
)
# Step 4: Subset the jittered points within Yangon window (Yangon owins)
quarter_conflict_ppp <- jittered_conflict_ppp[Yangon_Owin]
# Step 5: Rescale the ppp object to kilometers for KDE
jittered_ppp_km <- rescale(quarter_conflict_ppp, 1000, "km")
# Step 6: Calculate KDE using Gaussian kernel and bw.scott for bandwidth
kde_quarter <- density(jittered_ppp_km, sigma = bw.scott, edge = TRUE, kernel = "gaussian")
# Step 7: Store the KDE object for later plotting
kde_plot_list[[quarter]] <- kde_quarter
}Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points# Step 8: Plot all quarter KDE data using spatstat's plot function
par(mfrow = c(4, 4)) # Adjust rows and columns for multiple plots
for (quarter in names(kde_plot_list)) {
plot(kde_plot_list[[quarter]], main = paste("KDE for Yangon -", quarter))
}
5.5 Summary of Data sets
We have created a total of five datasets, each serving a distinct purpose in our spatial analysis:
ACLEDData_SF: This is the cleaned and formatted ACLED dataset, stored as an sf (simple feature) object. It contains the key event data, including attributes like location, event type, and time period, and will be used for further spatial analysis.M_SF: This dataset represents the Myanmar boundary as an sf object. It will be used as the base layer in the analysis, defining the spatial extent and providing a reference map for the region.Aggregated Level: This dataset is used for exploratory analysis and contains aggregated data for different event types, quarters, or other relevant categories. It allows for a broader overview before diving into more detailed spatial point pattern analyses.ACLED_ppp_corrected:This is the point pattern dataset created from the ACLED data. It includes the corrected coordinates (with jittering applied) and is stored as apppobject (planar point pattern) for use in spatstat. This object will be used for more in-depth spatial analyses, such as density estimation or clustering.myanmar_owin:This is the spatial window object representing the boundary of Myanmar. It defines the geographic limits for our spatial point pattern analysis and ensures that all event data is analyzed within this boundary.
5.6 getting the bandwidth
6.0 Deriving the Quarterly Kernel Density Estimation Layers
6.1 Setting the band widths
# Load necessary libraries
library(spatstat)
library(dplyr)
library(sf)
# Step 1: Define the function to create jittered ppp objects for conflict data
process_quarter_conflicts <- function(region_sf, region_window, region_name, data_sf) {
# Extract unique quarters
quarters <- unique(data_sf$quarter_year)
# Initialize empty lists to store ppp and KDE objects
region_quarters_conflict <- list()
kde_plot_list <- list()
# Step 2: Loop over each quarter and process conflict data
for (quarter in quarters) {
# Filter the data for the current quarter
quarter_data <- data_sf %>%
filter(quarter_year == quarter)
# Extract coordinates for this quarter
quarter_coords <- st_coordinates(quarter_data)
# Create a ppp object for this quarter
quarter_ppp <- ppp(
x = quarter_coords[, 1],
y = quarter_coords[, 2],
window = region_window
)
# Step 3: Remove rejected points that fall outside the window
valid_points <- inside.owin(quarter_ppp$x, quarter_ppp$y, region_window)
quarter_ppp <- quarter_ppp[valid_points]
# Step 4: Jitter duplicates
is_duplicate <- duplicated(quarter_ppp)
jittered_x <- quarter_ppp$x + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
jittered_y <- quarter_ppp$y + ifelse(is_duplicate, runif(npoints(quarter_ppp), -0.1, 0.1), 0)
# Create a jittered ppp object with the new coordinates
jittered_ppp <- ppp(
x = jittered_x,
y = jittered_y,
window = region_window
)
# Step 5: Store the jittered ppp object for later use
region_quarters_conflict[[quarter]] <- jittered_ppp
# Step 6: Rescale to kilometers and calculate KDE
jittered_ppp_km <- rescale(jittered_ppp, 1000, "km")
kde_quarter <- density(jittered_ppp_km, sigma = 100, edge = TRUE, kernel = "gaussian")
# Step 7: Store the KDE object
kde_plot_list[[quarter]] <- kde_quarter
}
# Return the processed ppp and KDE objects
return(list(
"ppp_list" = region_quarters_conflict,
"kde_list" = kde_plot_list
))
}
# Step 2: Define Yangon-specific variables
Yangon_ACLED_bbox <- st_bbox(Yangon_ACLED_Data_Sf)
Yangon_ACLED_window <- owin(
xrange = c(Yangon_ACLED_bbox["xmin"], Yangon_ACLED_bbox["xmax"]),
yrange = c(Yangon_ACLED_bbox["ymin"], Yangon_ACLED_bbox["ymax"])
)
# Step 3: Process the conflict data for Yangon
results <- process_quarter_conflicts(
region_sf = Yangon_Sf,
region_window = Yangon_Owin,
region_name = "Yangon",
data_sf = Yangon_ACLED_Data_Sf
)Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points
Warning: data contain duplicated points# Step 4: Plot the results
# Plotting the jittered conflict points for each quarter
par(mfrow = c(4, 4)) # Adjust rows and columns for multiple plots
for (quarter in names(results$ppp_list)) {
plot(results$ppp_list[[quarter]], main = paste("Yangon Conflict Events -", quarter))
}
# Plotting the Kernel Density Estimate (KDE) for each quarter
par(mfrow = c(4, 4)) # Adjust rows and columns for multiple plots
for (quarter in names(results$kde_list)) {
plot(results$kde_list[[quarter]], main = paste("KDE for Yangon -", quarter))
}
Mandalay_ACLED_bbox <- st_bbox(Mandalay_ACLED_Data_Sf)
Mandalay_ACLED_window <- owin(
xrange = c(Mandalay_ACLED_bbox["xmin"], Mandalay_ACLED_bbox["xmax"]),
yrange = c(Mandalay_ACLED_bbox["ymin"], Mandalay_ACLED_bbox["ymax"])
)
# Step 3: Process the conflict data for Mandalay
results <- process_quarter_conflicts(
region_sf = Mandalay_Sf,
region_window = Mandalay_Owin,
region_name = "Mandalay",
data_sf = Mandalay_ACLED_Data_Sf
)Warning: 139 points were rejected as lying outside the specified windowWarning: 123 points were rejected as lying outside the specified windowWarning: 152 points were rejected as lying outside the specified windowWarning: 157 points were rejected as lying outside the specified windowWarning: 120 points were rejected as lying outside the specified windowWarning: 139 points were rejected as lying outside the specified windowWarning: 124 points were rejected as lying outside the specified windowWarning: 113 points were rejected as lying outside the specified windowWarning: 149 points were rejected as lying outside the specified window
Warning: 149 points were rejected as lying outside the specified windowWarning: 234 points were rejected as lying outside the specified windowWarning: 161 points were rejected as lying outside the specified windowWarning: 210 points were rejected as lying outside the specified windowWarning: 51 points were rejected as lying outside the specified window# Step 4: Plot the results
# Plotting the jittered conflict points for each quarter
par(mfrow = c(4, 4)) # Adjust rows and columns for multiple plots
for (quarter in names(results$ppp_list)) {
plot(results$ppp_list[[quarter]], main = paste("Mandalay Conflict Events -", quarter))
}
# Plotting the Kernel Density Estimate (KDE) for each quarter
par(mfrow = c(4, 4)) # Adjust rows and columns for multiple plots
for (quarter in names(results$kde_list)) {
plot(results$kde_list[[quarter]], main = paste("KDE for Mandalay -", quarter))
}